Repurposing and computational design of PARP inhibitors as SARS-CoV-2 inhibitors

Coronavirus disease 2019 (COVID-19) is a recent pandemic that caused serious global emergency. To identify new and effective therapeutics, we employed a drug repurposing approach. The poly (ADP ribose) polymerase inhibitors were used for this purpose and were repurposed against the main protease (Mpro) target of severe acute respiratory syndrome Coronavirus 2 (SARS-CoV-2). The results from these studies were used to design compounds using the ‘Grow Scaffold’ modules available on Discovery Studio v2018. The three designed compounds, olaparib 1826 and olaparib 1885, and rucaparib 184 demonstrated better CDOCKER docking scores for Mpro than their parent compounds. Moreover, the compounds adhered to Lipinski’s rule of five and demonstrated a synthetic accessibility score of 3.55, 3.63, and 4.30 for olaparib 1826, olaparib 1885, and rucaparib 184, respectively. The short-range Coulombic and Lennard-Jones potentials also support the potential binding of the modified compounds to Mpro. Therefore, we propose these three compounds as novel SARS-CoV-2 inhibitors.


Strategy to design new compounds.
To design a novel compound, the target-ligand active site of Mpro was scrupulously examined, as it had a large volume. After careful visual analysis using DS, we identified two empty sites with olaparib (site 1 and site 2) and one empty site with rucaparib (site 1) and called them R 1 (Supplementary Fig. 2) .
At the identified R 1 position, we applied the "Grow Scaffold" approach, which has inbuilt libraries that facilitate the addition of groups to the small molecule ( Supplementary Fig. 3). The olaparib site 1 generated a total of 3,201 compounds, site 2 generated 304 compounds; and site 1 of rucaparib generated 437 compounds, accounting for 3942 compounds after these identified sites were subjected to Grow Scaffold modules (Fig. 2a). The simplified molecular-input line-entry system (SMILES) representations of these compounds are provided in (Supplementary Table 2).
The 3942 compounds were subjected to Lipinski's rule of five, enabling the 'Filter by Lipinski and Veber Rules' protocol. According to the "Rule of Five", a drug-like molecule should have maximum five hydrogen bond donors, 10 hydrogen bond acceptors, a molecular weight less than 500 Dalton, and a LogP value less than five. This process resulted in the production of 19 compounds. These compounds were subsequently docked to Mpro to delineate their binding affinities (Fig. 2b). www.nature.com/scientificreports/ Step 2: Binding affinity studies and clustering of small molecules. In this step, all 19 modified compounds were docked into the Mpro-binding pocket. The results revealed that three compounds displayed higher binding affinities than their parent structures (Table 1 and Supplementary Table 1). Manual clustering of the compounds was performed to understand the prospective binding modes of the small molecules (Fig. 3).
The pose with the best docking score from the largest cluster was further examined using molecular dynamics simulation (MDS). The two Olaparib-modified compounds were well aligned within the binding pocket of the target. Olaparib 1826 formed 27 poses in a cluster (Fig. 3a), whereas olaparib 1885 formed 18 poses in a cluster (Fig. 3b). Whereas, rucaparib184 formed a cluster with 7 poses (Fig. 3c). These results indicate that olaparib 1826 may have a higher affinity for Mpro, followed by olaparib 1885 and rucaparib184.  www.nature.com/scientificreports/ The three compounds demonstrated satisfactory results, with higher docking scores than the parent compounds (Table 1 and Supplementary Table 1). The parent compounds, olaparib and rucaparib, generated dock scores of 50.38 kcal/mol and 58.80 kcal/mol, respectively. The modified olaparib compounds, olaparib 1826 and olaparib 1885, demonstrated a dock score of 58.89 kcal/mol and 57.51 kcal/mol, respectively. Similarly, the modified rucaparib compound, rucaparib 184, had a dock score of 58.80 kcal/mol. These findings suggest that the modified compounds demonstrate a stronger affinity for Mpro than the parent structure (Table 1 and  Supplementary Table 1).
From the largest cluster, the compounds with the best docking scores were upgraded to molecular dynamics simulation (MDS) to understand the binding potential of the modified compounds to the binding pocket of the target. MDS analysis. MDS was performed for 100 ns to analyse the binding potential of the ligands at the binding pocket of the targets. The MDS analysis was based on the root mean square deviation (RMSD), radius of gyration (Rg), root mean square fluctuation (RMSF), and hydrogen bond number.
Stability analysis by RMSD. The protein backbone demonstrated stability 42 throughout the simulation. All the systems converged well and exhibited RMSD below 0.3 nm. The system with olaparib 1826 was largely stable during MDS evolution. At ~ 11,000 ps, a minute surge in the RMSD profiles was noted, and thereafter the system remained steady with an average of 0.16 nm (Fig. 4a). The RMSD profiles of olaparib 1885 and rucaparib 184 were stable throughout the MDS evolution without any noticeable variations, with an average of 0.17 nm and 0.2 nm, respectively (Fig. 4a).
Compactness analysis by Rg. Rg is "defined as the ratio of the protein accessible surface area to that of an ideal sphere of the same volume" 43 . This principle governs the protein-folding mechanism 43 . Our analysis revealed that Rg defines the distance between each atom of a protein and its centroid, thereby defining its compactness 44 . These three systems showed that the protein backbone was highly stable and compact. The readings exhibit a range of 2.20 nm to 2.28 nm with an average of 2.21 nm, 2.22 nm, and 2.22 nm for olaparib 1826, olaparib 1885, and rucaparib 184, respectively (Fig. 4b). This finding demonstrates that the systems were remarkably compact.
Fluctuation analysis by RMSF. The RMSF plots reveal fluctuations during the MDS run. All the systems were stable, with no major fluctuations. The systems were stable with RMSF below 0.3 nm (Fig. 4c). The average RMSF was 0.17 nm, 0.12 nm, and 0.13 nm for olaparib 1826, olaparib 1885, and rucaparib 184, respectively (Fig. 4c). This result indicates that each residue in the system was stable.
Number of hydrogen bonds. The number of hydrogen bonds between the protein and ligand was assessed during the entire simulation run. The presence of stable interactions suggested that ligands were held within the binding pocket throughout the simulation run. The three systems demonstrated hydrogen bond interactions during the entire simulation run, with an average of 2.0 for olaparib 1826, 1.3 for olaparib 1885, and 4.4 for rucaparib 184 (Fig. 4d). Rucaparib184 demonstrated a higher number of hydrogen bonds than the olaparib www.nature.com/scientificreports/ derivatives (Fig. 4d). This result implied that the compounds were firmly presented in the binding pocket of the protein.
Binding mode analysis. From the stable RMSD of the last 10 ns (~ 90,000-100,000 ps), representative complex structures were extracted and superimposed onto the X-ray crystal structure of Mpro. The results demonstrated that the ligands occupied a binding pocket similar to that of the co-crystallised ligand (Fig. 5), held by various intermolecular interactions. The modified compounds settled at different subsites of the target-binding pocket held by various interactions.
Key residue interaction of the modified compounds with the target residues. Olapar-ib1826. The modified compound olaparib 1826, formed hydrogen bonds with Gly143 and Gln189 (Fig. 9a). Interestingly, oxygen atoms were involved in these interactions. The distances between the atoms involved in hydrogen bonding were calculated. The Gly143: HN-O3 interaction was highly stable, with an average of 0.21 nm  www.nature.com/scientificreports/ throughout the simulation run (Fig. 6a). Similarly, the interaction between Gln189:HE21-O2 was also stable with an average of 0.28 nm, which indicates that the ligand tightly adheres to the binding pocket of the target (Fig. 6b). The key residue, His164, formed carbon hydrogen bonds to hold the ligand in the binding pocket. The residue Cys145 formed the π-alkyl interaction with the modified compound. The residues His41, His164, and Met165 adhered to the compound via alkyl interactions. His164 interacted with the fluorine atom of the ligand, thereby settling the ligand in the binding pocket of the target. Other residues, such as Met49, Phe140, Leu141, Asn142, Ser144, His163, Glu166, His172, Asp187, and Arg188 (Table 1) firmly held the compound in the binding pocket via van der Waals interactions (Fig. 9d).
Olaparib 1885. Olaparib 1885 formed hydrogen bonds with the Gly143, Ser144, and Cys145 residues (Fig. 9b). The hydrogen bond distance between Gly143: HN-O3 was stable, with no major aberrations. The average distance measured was found to be 0.2 nm (Fig. 7a). The hydrogen bond distance between Ser144: HN-O3 was analysed during the evolution of MDS (Fig. 7b). These findings reveal that the distance was largely stable, with a dip in the profile at 89,900 ps, which remained stable thereafter. Interestingly, the overall average distance was measured at 0.33 nm, while the last 10 ns were measured at 0.26 nm, suggesting that Ser144 strongly interacted with the ligand (Fig. 7b). Hydrogen bond interactions between Cys145: HN-O3 were stable throughout the simulation. There was a minor depression in the distance plot at 85,000 ps, which remained stable thereafter. While the overall distance was measured to be 0.31 nm, the average measurement for the last 10 ns was 0.23 nm (Fig. 7c). Generally, these three interactions are firm and hold the ligand in the binding pocket. Residues Asn142, His163, and His164 formed carbon-hydrogen bonds that firmly held the ligand at the active site. The key residues His41 and Cys145 have generated alkyl and π-alkyl interactions, thereby positioning the ligand at the binding site.
The hydrogen bond distances between the protein residue atoms and ligand atoms indicated that the interactions were stable during the simulation run. The Arg40: HE-O47 interaction has seen an elevation during the initial steps of the simulation run to 15,710 ps; however, it was soon stable thereafter with an overall average of 0.19 nm (Fig. 8a). Another interaction was observed between His41: HD1-O49. The distance profile displayed a stable interaction, although a few fluctuations were observed from 10 to 11 ns. The average distance was measured to be 0.25 nm during the evolution of the simulation (Fig. 8b). The hydrogen bond interaction distance between Gln192: HE21-N5 was also stable throughout the simulation. However, the average distance was projected to be 0.46 nm. This distance measured seems to be marginally weaker than the previous interactions that allowed the ligand to settle at the binding pocket (Fig. 8c). His164 and Glu189 residues adhered to the compound via carbon-hydrogen bonds. The key residues Met49, His41 and Met165 have prompted π-alkyl and π-π T shaped interactions with the ligand. The residue Cys44 formed an interaction with the fluorine atom, and Arg40 generated an attractive charge interaction, holding the ligand in the binding pocket of the compound. The residues Thr45, Ser46, Tyr54, Cys85, Glu166, Leu167, Pro168, Phe181, Thr190 and Ala191 (Table 1) have interacted with the ligand via the van der Waals interactions, thus positioning the ligand at the binding pocket (Fig. 9f). Analysis of the modified groups of the three compounds revealed that rucaparib 184 generated two hydrogen bonds with key residues.   www.nature.com/scientificreports/ than both olaparib 1826 and olaparib 1885. These findings suggest that the modified compounds formed thermodynamically firm bonds with Mpro, as reported earlier 45 .

Discussion
PARP proteins demonstrate structural similarity and function, with two riboses and two phosphates in a unit polymer 36 . PARP1 was discovered in 1963 36 . They perform critical functions including apoptosis, DNA damage response, and transcription modulation 46,47 . Their inhibitors exert anticancer activities 36,46,48 . In the current study, PARP inhibitors were used against COVID-19 target to specifically design putative inhibitors with affinity for SARS-CoV-2 Mpro. To identify effective therapeutics for COVID-19, PARP inhibitors were employed as starting structures. Here, we computationally designed compounds that demonstrate higher affinity for Mpro than their parent structures.
To design new molecules computationally, we first identified a target with a large active site volume. We then marked the unoccupied spaces within the binding pocket after molecular docking. Correspondingly, Mpro was identified as the target and PARP inhibitors were docked to it. Olaparib and rucaparib had higher docking scores for Mpro. Therefore, these compounds were used in subsequent experiments.
The Mpro-ligand complex was selected to design new compounds by identifying the empty space located in the binding pocket (unoccupied by the ligand). Accordingly, with olaparib, two sites were identified (site 1 and site 2) and one site was identified with rucaparib (site 1) (Fig. 11, Supplementary Fig. 2). These sites were marked R 1 and the Grow Scaffold module was initiated.
Molecular docking was initiated, followed by clustering analysis, to evaluate whether the new compounds induced conformational changes. The results revealed that the modified groups occupied positions at the desired sites, which aligned with the hypothesis. Furthermore, the binding affinity results showed that two compounds (olaparib 1826 and olaparib 1885) from olaparib site1 showed stronger affinity than site2 and occupied a defined position (Fig. 11), as did one compound from site3 (rucaparib 184). The olaparib-modified compounds showed better results against Mpro than the rucaparib-modified compounds. However, the modified compounds generated higher docking scores than the parent compounds.
Olaparib 1826 was obtained from the fragment library of organosilanes, reaction name yama coupling, and fragment name 2-(triethoxysilyl)ethylamine. Olaparib 1885 was obtained from the fragment library Grig-nardReagents, the reaction name Kumada Coupling, and the fragment name cyclobutylmagnesium chloride. Rucaparib 184 was synthesized using the fragment library acids with the fragment name monosodium L-aspartate dehydrate #2 and a reaction called esterification.
In the MDS studies, olaparib 1826 and olaparib 1885 maintained a similar binding mode to that seen in molecular docking, while rucaparib 184 changed dramatically from its initial position. However, all the compounds were accommodated at the active site of the target while adhering to key residue interactions (Fig. 12).
Molecular docking interactions revealed that with Mpro, olaparib 1826 had hydrogen bond interactions with the key residues Gly143, Ser144, Gln189, and Thr190. MDS retained hydrogen bonds with Gly143 and Gln189. Interactions with Gly143 have been reported previously 28,30,49 . Similarly, a hydrogen bond was also observed with Gln189 in the molecular docking and MDS pose, as reported previously [50][51][52] . Hydrogen bonds between these residues were also observed in the X-ray structures. In addition, several other interactions originating from the binding pocket of the target firmly hold the ligand at its active site. The distance measured for these residues www.nature.com/scientificreports/ was also within the acceptable length below 0.3 nm during the progression of MDS (Fig. 6). Moreover, targeting the catalytic dyad residues His41 and Cys145 is important for the developing of Mpro inhibitors 53 . Our results show that olaparib 1826 formed van der Waals interactions with His41, and alkyl interactions with Cys145. No change in the binding mode was observed in the molecular docking or MDS poses (Fig. 12a). The compound olaparib 1885 formed hydrogen bonds with Gly143, Ser144, and Cys145. Similar bonds have been observed in a previous study [54][55][56] . The residue Glu189 that formed a hydrogen bond with the docked pose, also formed van der Waals interactions after MDS. However, the remaining two hydrogen bonds with Gly143 and Ser144 were retained. Hydrogen bond interactions with Gly143 and Ser144 were also observed in the X-ray structure; and Cys145 formed van der Waals interactions. The binding modes of olaparib 1826 and olaparib 1885 did not significantly change (Fig. 12b).
Rucaparib 184 formed hydrogen bonds with Arg40, His41, and Gln192 after MDS which was different from the molecular docked pose. Among the three compounds, rucaparib underwent drastic changes from its initial structure and was buried deep in the active site (Figs. 11 and 12c). The interactions between Arg40 and His41 were stronger, with an acceptable distance measured below 0.3 nm. Various interactions with Gln192 have been reported earlier 57,58 . The key residue Met165 had formed π-alkyl interactions with olaparib 1885 and olaparib 1826 while with rucaparib it has formed an alkyl interaction as was seen in the X-ray structure. However, in the molecular docking pose, this residue generated a carbon-hydrogen bond (Supplementary Fig. 4). Several other key residues originating from different subsites 59,60 also prompted van der Waals interactions. In particular, Glu166 promoted van der Waals interactions with all the ligands (Table 1).
We also calculated the binding energies for the three complexes along with the co-crystallised compound. This was done by initiating the 'Calculate Binding Energies' tool available on the DS. This calculation allows for the assessment of the binding energy between the receptor and the ligand. The binding energy for the cocrystallised compound was estimated to be − 136.724 kcal/mol. The binding energy for olaparib 1826, olaparib 1885, and rucaparib 184 were calculated to be − 44.2554 kcal/mol, − 63.7637 kcal/mol, and − 86.8231 kcal/mol, respectively. The higher results for the co-crystallised ligand with respect to the molecular docking score and binding energy may be due to its larger size than the modified compounds.  www.nature.com/scientificreports/ The identified compounds also demonstrated an acceptable Lipinski's rule of five which was calculated using the Filter by Lipinski's and Veber. These compounds showed acceptable results ( Table 2). Since the compounds were novel, we estimated their ability to be synthesised by adapting the SWISS ADME 61 web tool and read it according to the synthetic accessibility score. The results ranged from 1(very easy) to 10 (very difficult) 61 . The synthetic accessibility score for olaparib 1826, olaparib 1885, and rucaparib 184 were 3.55, 3.63, and 4.30, respectively. Scores close to 5, indicated that these compounds may be easily synthesised ( Table 2). These elegant findings show the compounds olaparib 1826, olaparib 1885, and rucaparib 184 are possible Mpro inhibitors.
The novelty of the compounds was verified using SMILES as an input to PubChem 62 . These results indicate that these compounds were not yet used. We speculated that these designed compounds (Fig. 13) have not yet been synthesised and could be new compounds for the treatment of COVID-19.

Materials and methods
Selection of the ligands. PARP inhibitors were sketched using Biovia Draw v2017 and saved as a molfile (. Mol) format (Fig. 14) and were exported to the DS. The ligands were minimised using the "minimize ligands" protocol available with the DS. The CHARMM force field, was adapted using a Smart Minimizer algorithm. This was executed using 1000 steps of steepest descent with a RMS gradient tolerance of 3; subsequently, conjugate gradient minimisation was applied. www.nature.com/scientificreports/ Selection of the targets. The target selected for this study was Mpro (PDB ID:6LU7) from SARS-CoV-2 26 .
The protein was prepared by enabling the 'Prepare Protein' protocol available with the DS. This protocol prepares and checks the given protein and performs actions that include modelling the missing loop regions, standardising the atom names, dislodging the water molecules, inserting missing atoms in the incomplete residues 63,64 , deleting the alternate conformations, and protonating titratable residues using the predicted pKs. The active site for molecular docking was selected around the co-crystallised ligand for all atoms and residues at 13.82 Ǻ. This was done by enabling the tool 'Define and Edit Binding Site' that correspondingly creates a binding sphere. It has been reported that the active site contains subsites of origin of different residues 65 .
Prior to the initiation of molecular docking, the co-crystallised ligand N3 was dislodged and redocked into the selected binding pocket to ensure that the binding mode was reproduced. The results have shown that the ligand generated a similar binding mode as that of the co-crystallised ligand with an acceptable RMSD of 0.9 Ǻ ( Supplementary Fig. 5).
Binding affinity studies. Binding affinity studies were performed between Mpro and small molecules using the CDOCKER program available in the DS. Prior to the initiation of docking studies, the binding pockets of the targets were examined to determine the active sites. We aimed to identify a target with a large active site to facilitate the designing the SARS-CoV-2 specific inhibitors ( Supplementary Fig. 1).
Accordingly, Mpro was selected, and PARP inhibitors were docked into the active site. Correspondingly, the active site was chosen to be around the co-crystallised inhibitor N3. The selected ligands were used to generate 30 conformations. The best pose was selected after clustering the conformers to determine the best binding mode. From the largest cluster, the compound with the best docking score (binding affinity) revealing interactions with the key residues, was chosen.

Strategy to design the new compounds.
To design new compounds computationally, we examined the target ligand complex, after molecular docking, to locate the empty spaces. Subsequently, certain points are identified and marked as R 1. The "Grow Scaffold" module available with the DS was enabled which resulted in a total of 3942 compounds. These compounds were docked into Mpro to estimate their binding affinities after a drug-like assessment.
The 'Grow Scaffold' module permits the user to accomplish reaction-based ligand listing inside the protein's active site by performing lead optimisation. It begins with the positioning of the scaffold at the receptor-binding site. The user can select a position that can act as a reaction vector, followed by the reactions and reagents to be used. Furthermore, the ligand "novelty" is computed by sorting and ranking the ligands by number of chain assemblies, number of double and aromatic bonds; and N, S, O atom count.

Molecular dynamics simulation (MDS) analysis.
The protein-LIG complex structures obtained from molecular docking were escalated to MDS to understand their binding potency at the binding pocket of the targets. GROningen MAchine for Chemical Simulations (GROMACS) v2016.6 was used to study the MDS 66,67 . A CHARMM27 all-atom force field was utilized 68 . The topologies of the ligands were acquired from SwissParam 69 and the topol.top file was updated accordingly. Subsequently, the systems were solvated in a dodecahedral water www.nature.com/scientificreports/ box using the TIP3P water model. The system was neutralised using counterions. Energy minimisation was performed in 50,000 steps to remove bad contacts and clashes. After successful energy minimisation, the protein and ligand were coupled using the gmx make_ndx command to progress through the two-step equilibration process. The first step of equilibration was (constant number of particles, volume, and temperature) NVT equilibration, which was performed for 100 ps using a V-rescale thermostat at 300 K. The second equilibration was conducted with a constant number of particles, pressure, and temperature (NPT) for 100 ps using a Parrinello-Rahman barostat to monitor the pressure at 1 bar. During equilibration, the protein backbone was restrained and the non-protein was permitted to wobble. The long-range electrostatic interactions were evaluated by the Particle Mesh Ewald (PME) method and the short-range interactions and interactions by van der Waals were calculated after applying upper limit of 9 Å and 14 Å, correspondingly. MDS progressed under periodic boundary conditions for 100 ns. The corresponding analysis was conducted using various GROMACS tools and visual molecular dynamics (VMD) 42,70 .
Analysis of the trajectory. The different tools available to analyse the MDS results, and GROMACS were utilized 71 . RMSD was assessed using the gmx rms. These calculations provide knowledge of the deviation present, if any, in the protein from the initial to final conformation during the MDS run. Logically, the lower the RMSD, the greater the stability of the protein 72 . Rg, which determines the compactness of the protein, was studied using gmx gyrate 73 . The RMSF of the protein residues was determined using the gmx rmsf tool. Here, the fluctuations in each residue of the protein were examined. The RMSD, Rg and RMSF were computed for the protein backbone. Using the gmx hbond, the number of hydrogen bond interactions between protein and ligand atoms were evaluated during MDS progression. Furthermore, the interaction energy was calculated using the gmx energy, which computes both the short-ranged Coulombic and Lennard-Jones energy interactions 45 .

Conclusion
PARP inhibitors are well-known for their anticancer properties. Recently, in vitro studies also demonstrated their anti-COVID-19 properties. In the present study, we aimed to design new PARP inhibitors with high affinity for the COVID-19 Mpro. Correspondingly, two compounds originating from olaparib and one compound originating from rucaparib displayed high binding affinities for Mpro. This indicates their potential use as SARS-CoV-2 inhibitors. These compounds also demonstrated favourable synthetic accessibility scores. Additionally, our method could aid the scientific community to design new compounds and discover new horizons for drug discovery against SARS-CoV-2.

Data availability
The generated dataset resulted from the 'Grow Scaffold' are provided as the Supplementary Table 2. A SMILES representation of the compounds is provided.